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Abstract 

In mathematical finance and other applications of stochastic processes, it is frequently the case that 
the characteristic function may be known but explicit forms for density functions are not available. The 
simulation of any distribution is greatly facilitated by a knowledge of the quantile function, by which 
uniformly distributed samples may be converted to samples of the given distribution. This article analyzes 
the calculation of a quantile function direct from the characteristic function of a probability distribution, 
without explicit knowledge of the density. We form a non-linear integro-differential equation that despite 
its complexity admits an iterative solution for the power series of the quantile about the median. We give 
some examples including tail models and show how to generate C-code for examples. 



1 Introduction 

The construction of Monte Carlo samples from a distribution is facilitated if one has a knowledge of the 
quantile function w(u) of a distribution. If F(x) is the cumulative distribution function associated with a 
continuous density, then the quantile w(u) is the solution of the equation 

F(w(u))=u. (1) 

A knowledge of the function w(u) makes Monte Carlo simulation straightforward: given a random sample U 
from the uniform distribution, a sample from the target distribution characterized by f(x),F(x) is 

X = w(U) . (2) 

This approach is especially useful when the [/-variables arise as the components of a sampled copula, as then 
one needs the quantile functions of the marginals to create the marginal samples from the copula sample. But 
in general the method is useful, irrespective of whether copulae are involved. 

A full differential theory of quantiles based on a density function has been given recently by Steinbrecher 
and Shaw [TT]. The method is particularly easy to apply analytically when the density is of Pearson form 
and various power series solutions may be developed. This new paper addresses the case where the density is 
not known explicitly but where one just has its characteristic function. Such cases are of increasing interest, 
especially due to the role of jump diffusions and Levy processes in mathematical finance. The methods 
developed here may also be applied to the cases of known density of more complicated form, as working via 
the characteristic function can turn out to be straightforward. We shall be able to give a detailed development 
here for the case where all integer moments of the characteristic function exist, i.e. where the probability 
density function is analytic about the origin. 
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2 Quantile— characteristic function relationships 

If f(x) is the probability density function for a real random variable X, the first order quantile ODE is obtained 
by differentiating Eqn. (1), to obtain: 

dw(u) 

/(w(u)) ^T = 1 ' (3) 

where w(u) is the quantile function considered as a function of u, with < u < 1. Applying the product rule 
with a further differentiation we obtain: 

au z V ait 



+ />(«))( -7^1 -0- ( 4 ) 



This may be reorganized into the second order non- linear ODE given by Steinbrecher and Shaw |llj : 

d w(u) . ,-.-.{ dw(u) s 



where 

H(w) = -~log{f(w)}. (6) 
dw 

and the simple rational form of H(w) for many common distributions, particularly the Pearson family allows 
analytical series solutions to be developed [11]. Changes of variable may also be introduced, leading to a 
useful differential characterization of Cornish-Fisher expansions and candidates for quantile representations 
optimized for GPU computation [5]. 

Now suppose we do not have an explicit form for f(w), but rather only have a characteristic function <j)x{t) 
defined by 

4>x{t) = E[e ltx ] (7) 
So the theoretical density / and (j>x are related as Fourier transform pairs with the conventions 

/OO 1 poo 

f(x)e ltx dx , f(x) = — / 4>x(t)e- ltx dt (8) 

Combining these definitions allows us to write down the following three links between the quantile function 
and the characteristic function. First, the definition of the characteristic function can be written as 

/oo 
e ltx dF(x) (9) 
-oo 

where F(x) is the cumulative distribution function associated with f(x). We can then change variables to 
u = F(x), with x = w(u), to obtain a first quantile-characteristic link in the form: 

<j) X (t)= f e ltw ^du (10) 
Jo 

Second, by using the Fourier inversion theorem, we have, from Eqn. (|3]), 

dw(u) 



(t)x(t)e~ Uw( - u) dt = 2n, (11) 



du 

Third, and in the spirit of our previous approach, a further differentiation gives us 
d 2 w(u) f°° u ,-itwM^ Jdw(u)^ 2 



4> x (t)e- itw ^dt -if dv M) r t x (t)e- Uw ^dt = 0. (12) 
V du J J_ x 



du 2 

These last two relationships may be combined to give a fourth identity in the form: 
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This last equation will be the basis of our further analysis, though it is possible other efficient approaches 
may be developed from other such identities. We call it the characteristic quantile equation, or CQE. The 
CQE is a direct link between the characteristic and the quantile. If we can solve the CQE we can significantly 
short-cut what would be the normal three-stage route from the characteristic function to the quantile, and 
obtain Monte Carlo samples, as illustrated in the diagram below. 

, / i \ Inverse Fourier transform 



CQE solution? 

w(u) < 

Functional Inversion 

Sample uniform U j 

X = w{U) 

It is quite often the case the the first step is analytically intractable, and if it does succeed the third step of 
functional inversion may be difficult. The diagonal route from the characteristic function direct to the CDF 
F(x) has been considered, and we turn to this next. 



fix) 

^Integration 

F(x) 



3 From the characteristic function to the distribution function 

The diagonal route from <fcx(t) to F(x) is a well-trodden one. We refer to reader to Shephard (1991) [3] for 
a definitive survey on the results related to the "Gil-Pelaez inversion" formula [2] ■ We will discuss this result 
for two reasons. First we will need a special case for a result that anchors our quantile functions. Second, we 
wish to bring attention to the danger of using this result outside a carefully defined context - this comes to 
light if one considers an essentially complex analytic picture. 

The starting point is the elementary observation that the CDF F(x) can be regarded as the convolution 
of the density f(x) with the Heaviside step function 9(x) that is unity for x > and zero otherwise. 



F(x)= / f(y)dy= / 6(x-y)f(y)dy (14) 

J — CO J — CO 

We therefore wish to exploit the convolution theorem, the Fourier transform of F is given by 

F(t) = <Px(t)9(t) (15) 

Let us work out 9(t) carefully: 

/•CO /"CO 

0(f) = / 9(x)e ltx dx = / 6(x)e ltx dx = / e ltx dx (16) 



The last integral is convergent if and only if ^s(t) > 0, i.e. the transform exists and is holomorphic on the 
strict upper half-plane. In that case we obtain 

*" " FSj (17) 

In principle this means that an inversion integral involving 9(t) should be carried out along a horizontal 
contour in the upper half plane U = {z : ^s(z) > 0}. Now consider f(t). This might not in fact be holomorphic 
(we will return to this later), but if it is it will normally be holomorphic on a horizontal strip D = {z : B < 
3(z) < T}. Provided D D U is non-empty we pick a horizontal contour C € D PI U and the inversion and 
convolution theorems give us 

m - i / (is) 



2tt J c t 



Note that this whole idea requires that D n U is non-empty. Furthermore, any attempt to move C to the real 
axis might fail - this is not due to the presence of the factor 1/t, which is in fact easily managed. Two cases 
of notable interest in the wider context of option-pricing are 
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• /i(x) = (e x — 1) + : the reader may verify by elementary integration that D = {z : 1 < 3(2)} and the 
real axis is well below D n U = D; 

• = (1 — e x ) + : the reader may verify by elementary integration that D — {z : ^s{z) < 0}, then 

Dnu = 0. 

In these two examples f% = fa = l/(it— i 2 ) and the consideration of D is all-important. However, the reader is 
right to note that these examples are somewhat pathological with regard to convergence - neither corresponds 
to the transform of a genuine PDF, whatever their wider relevance. 

Now suppose that <j>x(t) is holomorphic in a neighbourhood of the real axis. Then D fl U contains C(e), 
the horizontal contour {z : z = x + is; —00 < x < 00} traversed from left to right. Then 

27r JC(e) t 

We now deform C(e) to the contour C(e) (the mousehole contour) that is the union of the intervals 

Ci = {t: -00 < t < e} 

C* 2 = {t : i = ee l6 ;ir > 6 > 0} (20) 
C 3 = {t : e < f < 00} 

and where Ci is traversed clockwise. Then we can write 

1 j 2^7 Cl (-»t) 2vry c2 (-it) 27r7 C3 (-it) 1 ' 

and using the calculus of residues on C2 , we have 

lim J- / ^ x{t)e ~ ltX dt = ^[-^(0)] - ^ x (0) = i (22) 

We now change variables in C3 with t —> —t and establish the result. 

Given that we have assumed <f> is holomorphic we can expand the numerator and observe that the limit e — > 
may be taken, finally giving us the Gil-Pelaez inversion formula: 

f(x) = - + — [°° (24) 

Note that this result is the "shadow" , evaluated along the real axis, of the result of equation (18) in a context 
where C may be pushed arbitrarily close to the real axis. It is simply not true that this can be done with a 
general Fourier transform, as previously noted. 

The category of Fourier transforms that we actually wish to deal with the set of transforms that arise 
from transforming probability density functions. This is a rather more subtle category, as the characteristic 
functions that arise are not necessarily in the category of functions holomorphic in a neighbourhood of the real 
axis. The contrast between the two-sided exponential distribution and the Cauchy distribution is useful. The 
two-sided exponential distribution has a transform with nice thick strip containing the real axis and bounded 
by two simple poles. The Cauchy distribution on the other hand has a transform of the form 

4>x(t) = (25) 

which is simply not holomorphic in a neighbourhood of the real axis (cf stable distributions in general). 
However, one can reconsider Eqn. (24) in the altogether different setting of transforms of actually probability 
density functions, where the transforms are considered along the real axis only. We will not discuss this matter 
further here, and refer the reader instead to the work by Shephard 9 , which establishes that in this context 
the G-P inversion still works. However, the temptation to use Eqn. (24) in any wider setting must be resisted, 
for the reasons already stated. 
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We now turn to the result that we shall need for the direct computation of quantile functions. We need 
this result to establish an origin for expansions. 

Theorem 1 - the zero quantile location 



The value of u satisfying 
is given by 



w(u Q ) = (26) 

2 27T J Q t 

This resolved is now easily proved by substitution of x = in the inversion formula. We note that as expected 
Mo = 1/2 when /, and hence 4>x is symmetric. 

3.1 Examples of using the zero location theorem 

Consider an elementary case where we have not exploited the location symmetry on the simple Gaussian case 
to set the mean to zero. Then for a distribution with mean /i we have 

cj) X {t) = e vt -* 2/2 (28) 
The location theorem gives us, following some first simplification: 

1 1 f°° f 2 /2 2sm< fit) , 1 1 J fi, \ Ar/ , 

where N{x) is the normal CDF. This of course is the correct result. 

3.2 Stable distribution 

Next, consider a stable distribution with 

<f>x(t) = exp{-|t| a (l - i/feign(t)*)} (30) 
and <I> = tan(a7r/2). The location theorem gives us 

1 1 f°° f c sm( /3<f>t a ) 

The integral evaluates to give us 



u = 1 - — arctan(/3$) . (32) 
2 Tree 

Note than when [3 — we obtain u = 1/2 and that when /3 = 1 we obtain u = 0. In the latter case the 
distribution only exists for non-negative X. When a < 1 zero is not an interior point of the support of X. 
Our location formula behaves in the right way and we can also see that as (i — > — 1 then u — > 1 . 

The point of these observations is that we can get a base point for a series expansion of the quantile by the 
computation of an integral. In the stable case we can see that the result is a simple "closed-form" formula. 



4 Formal solution of the characteristic quantile equation 



We have a non-linear integro-differential equation linking the quantile to the characteristic function. Inspection 
of the CQE indicates that we can differentiate it , then eliminate w"(u) from the right hand side by using 
Eqn. (13) again. Let's do this once: 



w 



"I P OO 1 /■ CO 

'(«) = K(u)] 3 — / it{-it)^ x {t)e- Uw ^dt + 2,[w'{u)] 2 w"{u)— / it^ x {t)e~ Uw{u) dt. (33) 

2?r J-oo 2lT J-oo 

Use of the CQE and some simplification leads to 
w '"(u) = [w'(uy 4 



^ J°° t^xity-^^dt + Zw'iu)-^^ 00 it4> x (t)e- itw ^dt 



(34) 
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Some experimentation lead us to write 

w (n) (u) = [w'(u)} n+1 P n [w'(u), w(u)} (35) 

where P n [w'(u), w(u)] is a polynomial in w'(u) of degree n — 2 whose coefficients involve integrals of w(u). 
The operation of differentiation and resubstitution gives the recurrence identity: 

d d 

P n+ i[x, w] = (n+ l)xP 2 [w]P n [x, w] + x 2 P 2 [w]—P n [x, w] + ^Pnix, w] (36) 



where 



1 r°° 

P 2 [x,w} = P 2 [w] = —J it<t> x {t)e- ltw dt (37) 
We can do an inductive verification of this relation and note its correctness when n = 2. 

4.1 The analytic assumption 

In carrying out the repeated differentiation of the CQE we are assuming that all derivatives of w(u) exist and 
that we may form a convergent power series with which to represent the solution. The operation of repeated 
differentiation on the quantile will call up integrals of the form 



1 f°° 

M k = — j t k mdt 



(38) 



and these are related to the derivatives of the density function (if they exist) by 

/< fc >(0) = {-i) k M k (39) 

from the Fourier inversion theorem. What our series solution is using the characteristic function to invert a 
Taylor series for the density about the origin without explicitly calculating it. In our initial solution to the 
problem, this series must exist in some form - other methods must be used if no such series exists. 

4.2 Symmetric distributions on the entire real line 

How might we use this to solve our problem? Let's consider an interesting special but rich case, where the 
density is symmetric about the origin, and the random variables X may take all real values. In the symmetric 
case if (1/2) = and we can consider a power series solution about the median. Inspection of the integrals 
obtained so far gives us w"(l/2) = and 



w'" 



1 r°° 

(1/2) = K(l/2)] 4 - / t^ x {t)dt (40) 



The normalization w'(l/2) is easily determined by the condition: 



1 



</(l/2) 2tt 



OO 

<t>x{t)dt (41) 



One observes that the formal solution for the quantile as a power series around the median is in principle 
completely determined by a knowledge of the even "moments" of the characteristic function , as these supply 
the coefficients of that power series. Furthermore they do so directly, without having to establish the density 
/ or CDF F. In this symmetric case we can write down the complete solution for the quantile as the series 

oo 2fc+l 

w(u) = v + w\l/2) ]T P 2k+ i K(l/2), 0] " , (42) 



k=l 



where v = w'(l/2)(u— 1/2). The main issue that remains is the computation of the P-cocfficicnts in terms of 
"moments" of the characteristic function, by the solution of the recurrence identity Eqn. (17) with the initial 
condition of Eqn. (18). One also expects such power series to become increasingly awkward as one moves into 
the tails, so we expect to have to supplement the model with a special treatment of the tail. 
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4.3 Locational invariance of the CQE 

The CQE has a notable, if rather obvious symmetry. It is invariant under the transformations 

4> x (t) -> e lta x (t) , w(u) ^w(u) + a. (43) 
This of course corresponds to a shift in the density 

/(*) - f(x - a) , (44) 

as one would expect. 

4.4 Scale invariance of the CQE 

The transformation 

t — > ct , (fix{t) — » 4>x{ct) , w(u) — > w(u)/c . (45) 

also generates a symmetry of the CQE. 

The scale and location invariance may often be used to transform the problem to one in standard form 
with fewer parameters. 

4.5 Asymmetric distributions on (—00,00) 

Once any relevant locational and scale transformations have been exploited, a distribution may remain asym- 
metric. In this case the choice of an origin for expansion must be resolved. In the symmetric case w(l/2) = 
is the median so we might consider whether the try to still expand about the median or about the point uo 
chosen so that w(u ) = 0. The advantage of the latter approach is that we can do the expansion again in 
terms of simple "moments" of the characteristic function, and we can compute an expression for u . We have 
the zero quantile location as a special case of the Gil-Pelaez inversion formula 

4.6 Formal solution for the asymmetric case 

Having found u using the methods just described, we then have, as before, 



W 



k 

{u) = v + w'{u )Y J Pk[w'{u )^] V — (48) 



k=2 



where v = w'(u )(u — u ). 



4.7 Exponential asymmetry 

For a significant class of distributions asymmetry in the system is introduced by an exponential scaling. That 
is, for a real parameter j3, 

f{x-,P) = C ^ ) f s { X y* (49) 
where c() is a normalization function and fs is a symmetric case. It is evident by the shift theorem that 

m0) = ^Mt-0) (50) 
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and hence the quantile must be determined only by integrals of the symmetric characteristic function ips- We 
can make this explicit by a complex shift of the integration contouiQ in the CQE with the change of variables 
t = p + if}. We obtain the modified CQE linking the full quantile with the symmetric ips'- 

5 Series solution of the CQE - the symmetric case 

Let us summarize the problem in the symmetric case. Given a characteristic function <f>x(t), we need to solve 
the iteration scheme: 

d d 

P n+1 [x, w] = (n + i)xP 2 [vj]P n [x, w] + x 2 P 2 [w] — P n [x, w] + ^P n [x, w] (52) 

ox Ow 

with the initial condition: 

1 r°° 

P 2 [x,w] = P 2 [w] = — / it<t> x {t)e- ltw dt (53) 



2tt . 

The power series of the quantile function about the median w{\/2) — is then 



» y 2k+l 



w(u) = v + w'(l/2) > P 2k+1 [ w '(l/2),0}— — (54) 



£ — -(2* + 1)1 



where v = w'(l/2)(u — 1/2) and 

1 1 



.'(1/2) -2.' (55) 

The solution of the iteration is a rather mindless process that is best automated in a symbolic computation 
environment, and we employed Mathematica to carry out the analysis. We summarize the results. Let 

1 f°° 

E k =^- t 2k cf x (t)dt (56) 
2tt J^oo 

denote the fc'th normalized even moment. We have already observed that 

P 3 [x,0]=E 1 (57) 

Symbolic iteration then yields further terms as follows, where we abbreviate Pk[x, 0] = pk, 

p 5 = 10xE 2 - E 2 

p 7 = 280x 2 E 3 - 56xE 2 E 1 + E 3 

p 9 = lU00x 3 Ef - A620x 2 E 2 E 2 + x (l26£f + 120E 1 E 3 ) - E± 

p n = UQUQQx A El - 560560a; 3 E 2 E\ + x 2 (l7160E 3 E 2 + 36036£f £i) - x (792E 2 E 3 + 220E 1 E i ) + E 5 
pis = 1 90590400a; 5 ^ - 95295200a; 4 £:2 J Ei + x 3 (3203200£ 3 £i + 10090080£f £i) 

- x 2 (126126S| + 360360^1 E 3 E 2 + 50050E?£ 4 ) + x (1716-Bf + 2002E 2 E i + 364 J B 1 ^ 5 ) - E 6 
pi 5 = 36212176000a; 6 i; i 7 - 21727305600a; 5 £; 2 £; 5 + x 4 (775975200^3^ + 3259095840^£; 3 ) 

+ x 3 (-13613600-B4-B 3 - 147026880£ 2 i;3£ 2 - 102918816£f £i) 

+ x 2 (l23760E 5 E 2 + 1166880^ E x + 1361360^ 2 ^4^i + 2450448^1 £ 3 ) 

+ x (-1UA0E 3 E 4 - 4368E 2 E 5 - 560^i^ 6 ) + E 7 

The display here of further terms would become rather unwieldy. In any case, such expressions are best stored 
symbolically for subsequent simplification for a particular distribution. We do note that this is a one-off 
computation that once done to high order in a symbolic computation environment can then be transferred to 
another computer environment for implementation along with computation of the characteristic moments. 



(58) 



1 We can do this by Cauchy's theorem, given that /3 must be constrained so that / remains L 1 integrable and the transform 
does not develop singularities 
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5.1 Testing the approach with the normal distribution 

This is one where a complete algebraic characterization of the median power series has been given . 

4> x (t) = e-* 2 / 2 (59) 

All the relevant integrals exist: 



E k = — (°° t 2n e' t2 ' 2 dt = -2 n -^T (n+- 
2n J„ oc 7r \ 2 



(60) 



and in particular 

to' (1/2) = V2tt (61) 

Then the use of the series above, Eqns. (27)-(31), and some simplification yields the normal quantilc in the 
form: 

^ , 1 fiL*/*L , t 127^(«-|) 7 , 4369^ (u-i) 9 



2w u-~ +-V27T 3 / 2 



u 



2/3 V 2 J 15%/2 315^2 11340^2 



// 1 \ 13\ (62) 



348077r n /2 ( u _ i) 11 // i 

1 , Q I I 



u — 



89100^2 2 

Many more terms can be computed using some computer algebra. The generation of this particular series is 
however best handled with the methods of where a purely algebraic recursion may be used instead of 
Eqn. (27). 

5.2 Testing on the Student distribution 

For a Student t distribution with n degrees of freedom the characteristic function is given by 

2 1 -fn^-5 |iK 2 i^ (y/n\t\) 
^ (t) = r ' ( ' §) 2W (63) 

See e.g. the detailed discussion by Hurst [5] for an elegant derivation of this. The normalization is 

- = -'(V2) = ^i^ 1 ■ (64) 
All the relevant characteristic moments exist in the form 

4 *n-^r(fc + i)r(fc + f + \) 

Ek = (65) 

The method leads to the Student quantile series as given in [7] and as extended in [TT] . 

6 Symmetric Distributions of interest 

We may consider distributions whose density is known, or not known, in explicit form. 
6.1 Symmetric stable distribution 

These have been extensively discussed - see Nolan. Employing location and scale invariance allows us to focus 
attention on the symmetric case in standard form: 

<Px{t) = (66) 

The tail behaviour is given to leading order by 

asin(7ra/2)r(q) 
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All the relevant integrals exist: 

E k = l°° t 2k e-^dt=^r(^^\ (68) 



and in particular 



*"'< 1/2) " rjiTT/^J < 69 » 



6.2 Symmetric Generalized Hyperbolic 

This is a case, SGH, where the density is known in closed form. Nevertheless, while repeated differentiation of 
the density is possible for certain parameters, it turns out to be much more tractable to work in momentum 
space. The SGH characteristic function is given by 

and is associated with the density function 



M ' V2nK x (a6) (V8 2 + x 2 /a) 1 / 2 - x 1 ' 



We assume that a > and 8 > 0. There are several special cases: 

• Student: A = —v/2, 8 = y/i/, a — > 0+. 

• Ordinary hyperbolic: A = 1; 

• NIG: A = -1/2; 

• Symmetric Variance Gamma: 8 — > + . 
The value of the density at the origin is given by 



f(0) = K *-lM aS ) = /Z^l/lM (70) 

J[ 1 V2nK x (a6) (<J/a)V2-A \J 2 nS K x (a5) [ 1 



(73) 



and so 

[2^5 K x (aS) 

W(0H V^^R) 

The characteristic moments are given by (exploiting the symmetry) 



F 1 r : , k ( ^ \ X/2 K x (SVa^) 



Making the change of variable p = \J a 2 + z 2 

P 1 f°° 1 2 2\(fc-l/2) A 

Ek= nJ a {P " a) ap K^S) (75) 

This integral can be evaluated by Mathematica |12j . which returns the value 

Ei ^r,* +V2 )(?r' 2 teM ™ 

subject to computer-generated constraints k > 0, a > 0, 8 > 0, A < 2. 
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6.3 A Harder Example: the Levy Stochastic Area 

The distribution of the Levy Stochastic Area (LSA) is of considerable interest from the point of view of 
stochastic analysis and high order Monte Carlo simulation. We base our analysis on the approach of Schmitz 
|10j . For a pair of Brownian motions spanning a time At the LSA is given by 

L(At) = / (W 1 (t)dW 2 (t) - W 2 (t)dW 1 (t)) (77) 



o 



Its characteristic function, conditional on a known value of R = W 2 (At) + W 2 {At), is known and given by 

4>l{z) = (/)x(z)<Py(z) 

a ( \ zAt 

<t>x{z) 



4>y(z) = exp 



sinh(zAt) (78) 
R 2 



2At 



[zAtcoth(zAt) - 1] 



where Wj(0) = 0. There are two components. When the path goes nowhere (loop) R = and this is given 
by the random variable X. Then there is a further contribution from Y. This split is convenient as X has a 
trivial quantile function: 



At 



ii 



Qx(u) = — log (79) 

7T V t — U / 

Readers should consult [TO] for details and an extensive discussion of this entity and further references. The 
difficult part is finding Qy(u). The characteristic function is clearly symmetric. The properties of the system 
are best understood by first extracting the time-scaling. From the definition L and X are both O(At) times 
some time-scale-independent random variable. So we set Y — AtP and R 2 = r 2 At, so that r is the distance 
gone by a Brownian motion in unit time. Making the change of variables shows that 

*<"> = s'(s) (80 » 



and the characteristic function of P is just 



b P (s) = cxp[-y(scoths - 1)] (81) 



and the associated density function is 



1 [°° r 2 
9p(p) = ^ exp[isp - — (scoths - l)]ds (82) 



— CO 



When s is small or p is large the hyperbolic function may be expanded to give 



<M S )~exp( g-J (83) 



so that samples of P may be approximately constructed with 

P-r^Z (84) 
where Z is a Gaussian random variable obtained by a quantile or other means. 

7 Implementation 

There are a number of issues to address when it comes to implementation. These include 
• The management of a symbolic differential recursion; 
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• The choice of a computer language in which to solve the problem; 

• Tail management; 

• Error analysis; 

• Implementation in legacy or specialized computer environments. 

The implementation of the differential recursion summarized in Section 4 appears to necessitate a modern 
symbolic computation environment and we have implemented it in a few lines of Mathematica |12j . But it is 
important to point out that this needs only to be done once and we have solved it for the symmetric case, 



the first few terms being given explicitly in Eqn. (58 1 up to k = 60? The computation becomes labourious for 
k > 50. The solution can be saved and then exploited. We can carry out the remain parts of the calculation 
entirely in Mathematica, as it can evaluate and store the characteristic moments and generate the relevant 
power series for the central portion of the quantilc. 

However, while this constitutes a complete solution for a truncated central power series, many potential 
users may need the ability to migrate the model to another computer language. For example, financial 
specialists may require C/C++ implementations. Statisticians might want an implementation in R. Some 
scientific and engineering applications will require a FORTRANXA, where probably XX > 90. We can 
consider trying to migrate the solution to another language at various stages. In this first implementation 
we will consider a late migration model, with C/CH — h as a target. This offers the compromise of getting 
Mathematica to do all the hard symbolic work, but we can output very simple C code to embed in a function 
to do the work in another program, once the parameters have been fixed. Migration one stage earlier would 
require calling a number of special functions in C, in particular various types of Bessel function. There are 
libraries to do this. Migration a further step back would require the implementation in C of the solution for the 
Pk- In principle this could be carried out with the CForm command used below in the late migration approach. 
We have not considered how the symbolic recursion itself might be solved directly in other languages. 

The error analysis is difficult as we do not know the benchmark, apart from special cases. We can give 
estimates based on this. Furthermore, if some form of expression for the cumulative distribution function 
(CDF), F, is available, we can certainly estimate the round-trip error 

RTE(u) = F[(Q(u)} - u (85) 

and from this, given a density /, we can estimate the quantilc (EQE) error as 

EQE{u) = {F[Q{u)\ - u)/f[Q{u)] (86) 

which is of course the Newton-Raphson correction. 

In each case the central series will almost almost require a tail model. The power series will work well in 
a region e < u < 1 — e but the region 1 — e > u > 1 and its mirror will need separate management. 

7.1 Solution of the differential recursion 

This takes place as follows. First we define P-i. 

P[2, x_, w_] := 1/(2 Pi) Integrated t Phi [t] Exp[-I t w] , {t, -Infinity, Infinity}] 

Subsequent terms are then, for the general case without symmetry, given by 
P[n_, x_, w_] := 

P[n, x, w] = Expand [n*x*P [2, x, w] *P [n - 1, x, w] + 

x~2 D[P[n - 1, x, w] , x] P[2, x, w] + D [P [n - 1, x, w] , w] ] 

This is the key calculation. Two other operations rewrite the integrals that result in terms of the Ej~ and apply 
symmetries. The program that implements all this is given in Appendix A, where we exploit the symmetry 
to give a program to just work out every other term. The results are stored in a file and we have computed 
as far as pn in symbolic form. 
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7.2 Example 1: Symmetric Stable, central series for quantile 



The symbolic computation of the pk may be combined with the evaluation of the and u/(0) given in Section 
6.1, to give the central power series for the quantile. The first few terms are given by: 



w(u) 



3 r(f) (« 



er i 



r(^)r(f)-ior(£) 2 ))(,-i) ; 
i2orm 7 




This expression is useful for some basic checking in the central zone but is not accurate enough for serious 
Monte Carlo simulation. One interesting check is to ask the computer to symbolically invert the series, and 
then differentiate the result. We obtain 



/(*) = 



2(?ra) 



247ra 



+ O (x 6 ) 



(88) 



Readers may recognize this as Bergstrom's central series [T] for the density function, which is a nice point of 
verification^] For high precision work, we need more terms. The stored internal representation in Mathematica 
may be used, and in particular to test the precision of the result for the known Gaussian and Cauchy cases. 
In the Gaussian case the standard normal quantile is l/y2 times the stable quantile with a = 2. We may 
compare the results with Mathematica' s internal high-precision routines. We will plot the precision in the form 



log 
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Power Series from Stable Distribution :a = 2 
\/2 x {Internal High — Precision Quantile) 



(89) 




Figure 1: Precision of the Gaussian stable quantile: 0.5 <u< 0.95 
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Figure 2: Precision of the Cauchy stable quantile: 0.5 < u < 0.95 



We note that the error is about machine precision level for 0.5 < u < 0.84 and then rises into the tail, 
remaining at an acceptable level until some point around 0.94. Similar results apply to the Cauchy case, where 



2 The use of Mathematical s InverseSeries command to compute the quantile from a density series may also be invoked, but 
is not practical for very large series. 
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Figure 3: Full precision of the Gaussian stable quantile: 0.5 < u < 0.95 

no factor of l/\/2 is needed, and the error grows slightly faster than in the Gaussian case. Note also that 
the machine precision oscillations in the plot are due to restricting attention to compiled results within the 
plotting routine - the series becomes much more accurate than one part in 10 16 as one approaches the median. 
We can illustrate this by asking for one of the comparisons to be re-done in arbitrary precision mode. The 
result is shown in Figure 3, with quad- and double-precision limits shown as horizontal lines at —32 and —16. 

7.3 The stable quantile with a = 3/2: migration to C/CH — h 

Having tested the symbolic stable quantile on known cases a — 2,1, we now consider the intermediate case 
a = 3/2. In this case we may take our symbolic representation and convert it into working C/C++ code. There 
are some helpful intermediate steps. First, we introduce variables v = 2u — 1 and w — v 2 . Second, we use the 
symbolic computation engine to write the system in the standard nested multiplication form for polynomials 
(truncated power series). In Mathematica this involves the use of the HornerForm function. Finally we use 
CForm to output the code. The output of this takes the form: 

v* (1.74002161967547716294123 + 

w* (0.648419685395586217984681 + 
w* (0.452196867009616298571941 + 
w* (0.371651133863068554240291 + 
w* (0.32860784309392901699825 + 

w* (0.302354425634672539752731 + 
w* (0.285017852611712836182585 + 
w* (0.272940090886219734608815 + 



w* 


'0 


264185174645675414204326 


+ 


w* 


'0 


257629942234930169112297 


+ 


w* 


'0 


25257672730254977798406 + 


w* 


co 


248568819726342542363224 


+ 


w* 


io 


245294721047639736371026 


+ 


w* 


:o 


242535016745509439912893 


+ 


w* 


io 


240131243639121611881269 


+ 


w* 




237966799943489229416831 


+ 


w* 


:o 


235954779696462330935446 


+ 


w* 


:o 


234029956041290960707096 


+ 


w* 


:o 


232143339065715759679999 


+ 


w* 


:o 


230258380102825285626612 


+ 


w* 


:o 


228348256742232151582891 


+ 


w* 


io 


226393883462430267805038 


+ 


w* 


:o 


224382419242756015826282 


+ 


w* 




2223061215967177252663 + 




w* 




220161445928088782138854 


+ 


w* 




217948321166808082853286 


+ 


w* 


[0 


215669553853982023082385 


+ 
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w* (0.213330327149878379304031 + 

w* (0.210937771052290678843475 + 

w* (0.208500586950509676634869 + 

w* (0.206028714471339648851334 + 

w* (0.203533032028819495783199 + 

w* (0.201025084985861029012365 + 

w* (0.198516837154295178958399 + 

(0.196020442690960745941278 + 

. 19354803642357645356947*w) *w) ))))))))))) ) 

))))))))))))))))))))) 



and gives a directly usable implementation in C/C++ that can be pasted into a simulation code. Full details 
of how this was generated are given in Appendix B and are also available on-line. 

7.4 Tail model for the stable quantile 

The power series computed thus far is a power series about the median, and as shown in the example precision 
plot for the Gaussian, loses precision in the tail. In each case we need to augment the quantile model by 
producing a separate model for the tail. We can of course add more and more terms to the power series, 
but ultimately we are confronted by the fact that usually such series should diverge as \2u — 1| — > 1 and a 
polynomial truncation of the series is not enough. 

There is an interesting mathematical question here. Should one try to develop a general theory based on 
some asymptotic analysis of the characteristic quantile equation? Or should one employ known information 
about particular cases? In this note we will take the latter approach, as there is a great deal of useful 
information available. 

In the symmetric stable case the relevant asymptotic series for the distribution function were established 
over 50 years ago by Bergstrom [T], based on earlier work on Laplace transforms by Pollard and others [?]. 
An easily accessible paper jj] has a formula (Eqn 2.9 of [I]) whose integral gives, for x — > oo, and a ^ 2, and 
after some simplification 

k=l v 7 

We have found that the formal inversion of this series with N = 4 provides a satisfactory tail model with a 
reasonably kinkless join to the power series. The result takes the form, for a < 2, 

w(u) ~ | — + c + ci(l - u) + ca(l - u) 2 \ (91) 

where the coefficients are given by 
_ r(a)sin(^) 

C— 1 

7T 

cos (f) T(2a) 

r(a) 

7TCSC 2 

Cl = - 



ttcsc 2 (^) (2r(a)r(3a)sin(2fi) - 3 esc (f) T(2a) 2 sin 2 (7ra)) 



(>2 



it 2 cot (f) esc (f) (6(cos(7ra) + l)r(2a) 3 — 3(2cos(7ra) + 1)T (a)T (3a)T (2a) + cos(7ra)r(a) 2 r(4a)) 

6r(a) 5 

(92) 



The testing of such representations for intermediate a requires comparison with specialist models for the stable 
distribution. In this case we made a comparison of the results with those published on the web by J. Nolan, 
at 

\protect\vrule widthOpt\protect\href {http : // academic2 . american . edu\string~ jpnolan/ stable/quantile . dat} 
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We found that for 0.1 < u < 0.9 the relative error between our results and Nolan's easily is less than 10~ 6 , 
giving high confidence in the power series. With the simple tail model the precision over the entire range is 
as shown in the figure. 

7.5 Non-analytic density functions 

In the detailed analysis given thus far, we have considered the analytic case. By no means all densities of 
interest fall into this category, and we must also treat cases where not all of the characteristic moments exist. 
For example, consider the (symmetric) variance gamma (VG) case. Taking the limit 5 — *■ + in the generalized 
hyperbolic model gives us the characteristic function 

*«> = (^) A 

and the density 

/ n>\ A+l/2 i 

/w=y 7^m a - i/2 w«m) 

The full VG model falls into the category of exponential asymmetry, with characteristic function 

*®={cP + {t-iw) 

A careful expansion of the density reveals that it comprises one power series in x 2 and a second power series 
in x 2 times a;' 2 * -1 ), so that a generalization of the methods developed here is needed for general A. There 
are logarithmic contributions when A = 1/2. In general rp ~ t~ 2X as t — > oo so we can also see that not 
all the characteristic moments exist. One approach to VG is to exploit a non-uniform base distribution, as 
discussed in [8. - in particular an exponential base is a convenient starting point for VG, hyperbolic and normal 
distributions. 



(93) 



(94) 



(95) 



8 Conclusions 

We have shown in principle how to establish a power series about the median for quantile functions char- 
acterized by a characteristic function linked to a smooth but possibly non-explicit density. This has been 
elucidated in sufficient detail for symmetric distributions based on a detailed symbolic solution of a non-linear 
integro-differential equation. Further work is in progress to treat asymmetric cases and characteristic functions 
corresponding to non-smooth densities, as well as tail representations. 
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Appendix A: Mathematica code for the symbolic recursion 

The actual program used for the symmetric case was as follows. Here we work out every other term, but start 
the recursion carefully. First we give the starting values: 

PP[2, x_, w_] := 
1/(2 Pi) Integrate [ 

I t \ [CapitalPhi] [t] Exp[-I t w] , {t , -Infinity, Infinity}] 

PP[3, x, w] = 

Expand [3*x*PP [2, x, w]*PP[2, x, w] + 

x~2 D[PP[2, x, w] , x] PP[2, x, w] + D[PP[2, x, w] , w] ] ; 



Next we give a second order iteration based on a repeated application of the first order form: 

PP[n_, x_, w_] := 
PP[n, x, w] = 
Expand [n*x* 
PP[2, x, 

w]*((n - l)*x*PP[2, x, w]*PP[n - 2, x, w] + 
x~2 D[PP[n - 2, x, w] , x] PP[2, x, w] + 
D [PP [n - 2 , x , w] , w] ) + 
x~2 D[((n - l)*x*PP[2, x, w] *PP [n - 2, x, w] + 
x~2 D[PP[n - 2, x, w] , x] PP[2, x, w] + 
D [PP [n - 2 , x , w] , w] ) , x] PP [2 , x , w] + 
D[((n - l)*x*PP[2, x, w]*PP[n - 2, x, w] + 
x~2 D[PP[n - 2, x, w] , x] PP[2, x, w] + 
D[PP[n - 2, x, w] , w]) , w]] 

Next we give the rules that set to zero terms that must vanish by symmetry: 

rulesa = Table [ 
Integrate [ 

I*t~(4 k + 1) *\ [CapitalPhi] [t] , {t , -Infinity, Infinity}] -> 
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0, -Ck, 0, 100}] ; 
rulesb = Table [ 

Integrate [-I*t~ (4 k - 1) *\ [CapitalPhi] [t] , {t , -Infinity, 
Infinity}] -> 0, {k, 1, 100}]; 
rules = Join[rulesa, rulesb] ; 

Then we give some results to write integrals in terms of the 

srules = Join[ 

Table [Integrate [ 

t~(2 k) \ [CapitalPhi] [t] , {t , -Infinity, Infinity}] -> 
2 Pi EE [k] , {k, 1, 200}] , 
Table [Integrate [-t~ (2 k) \ [CapitalPhi] [t] , {t , -Infinity, 
Infinity}] -> -2 Pi EE [k] , {k, 2, 200}]]; 

We initialize an array for solution: 

parr = Array [p, 200] ; 

We run the iteration as far as time allows, writing a file every ten steps, or every iteration after 60, to secure 
progress made so far. 

Do[dummya = PP[k, x, w] /. w -> 0; 
dummyb = (dummya /. rules) /. srules; p [k] = Collect [dummyb , x] ; 
If [k >= 7, (PP[k - 2, x, w] = 0)] ; 

Print [k, " Memory in Use = ", N [MemorylnUse [] /10~6] ] ; 
If [Mod [k + 1, 10] == I I k > 60 , 
Save[ToString[k] <> "pcoeff", parr]], {k, 3, 101, 2}] 

In our case we ran the code on a 2.8GHz Mac Pro and secured terms up to p^. We would be interested to hear 
about more efficient solutions. However, this many terms is sufficient for some serious practical applications, 
and we worked with the file 71pcoef f for the computed examples. 

The Mathematica notebook CharacteristicQuantileSeriesGen to generate the coefficients and the out- 
put files xxpcoef f s are available from the links at the web site 

www.mth. kcl . ac.uk/~shaww/web_page/papers/charquantiles/ 

Appendix B: Generation of C/CH — |- code for the stable quantile 

Here we show how to generate C/C++ code in the late migration approach for the case of the stable quantile. 
The Mathematica notebook CharacteristicQuantileExamples to do this is available from the web site 

www.mth. kcl . ac . uk/~shaww/web_page/papers/charquantiles/ 

The operations needed are as follows. First we load the general scries data: 
« "71pcoeff" 

Next we define the characteristic moments, w'(0) and some rules: 

g[k_] = 2 Integrate [t" (2 k) *Exp [-Abs [t] ~\ [Alpha] ] , {t , 0, Infinity}, 

Assumptions -> {\ [Alpha] > 0, k > -1/2}]; 
stablerules = Table[EE[k] -> 1/(2 Pi) g[k] , {k, 1, 60}]; 
wdash = l/( 1/(2 Pi) g[0]); 

The definitions of the rules rulesa, rulesb, rules, srules given in Appendix A are also re-loaded, and 
the series is then constructed in two steps: 

series = Table [(( (parr [ [2 m + 1]] x/(2 m + 1)! ) /. rules) /. 

srules /. stablerules), {m, 1, 35}]; 
StableSeries = series /. x -> wdash; 
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Creating the Gaussian example 

This is now easily done: 

GaussStableSeries = (series /. x -> wdash /. \ [Alpha] -> 2); 
wdashGauss = (wdash /. \ [Alpha] -> 2); 
GaussStableQuantile [u_] = 

1/Sqrt[2] (wdashGauss* (u - 1/2) + 

Sum [GaussStableSeries [ [k] ] * (wdashGauss* (u - 1/2)) "(2 k + 1) , {k, 
1, 35}]); 

Creating the a = 3/2 example 

This is now easily done: 

threehalf StableSeries = (series /. x -> wdash /. \ [Alpha] -> 3/2); 
wdashthreehalf = (wdash / . \ [Alpha] -> 3/2) ; 
ThreeHalf StableQuantile [ 

u_] = (wdashthreehalf *(u - 1/2) + 
Sum [threehalf StableSeries [ [ 

k]]* (wdashthreehalf* (u - 1/2)) "(2 k + 1) , {k, 1, 35}]); 

The function ThreeHalf StableQuantile may then be used within Mathematica. Within the late migration 
approach, we convert to C/C++ by first creating a suitable form for export. The following code accomplishes 
three things 

• works in terms of 2u — 1 in order to get radius of convergence unity and better control of the scale of 
the coefficients; 

• uses variables v = 2u — 1, w = v 2 ; 

• converts to nested multiplication form. 

export ThreeHalf = 
HornerForm [ 

N [ (wdashthreehalf *v/2 + 

Sum [threehalf StableSeries [ [ 

k] ]* (wdashthreehalf* (v/2))~ (2 k + 1) , {k, 1, 35}]), 24]] /. 
v~2 -> w; 

Then to generate the C/C++ code listing in Section 7.3 , we merely apply 
CForm [exportThreeHalf ] 
to obtain the results shown. 



